home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11tl.lha / lbl / hips / sources / 3d-tools / 3dog.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-06  |  25.5 KB  |  941 lines
  1. /*    Copyright (c)    1990    Jin Guojun
  2. %
  3. % 3dog.c:
  4. %    filter an image by applying Difference of Gaussians mask.
  5. %    The input is in byte or float format, and the output is in
  6. %    floating point or integer format.
  7. */
  8.  
  9. char    usage[]="\n\
  10. 3dog [-A esigma [masksize [ratio ]]] [-p #] [-i [-c]] [-m] [-w]\n\
  11.     [-b #] [-n #] [-s #] [-g #] [-flvt] [< input_seq] [> out_seq]\n\
  12. \n\
  13. esigma [1.0] is the standard deviation of the \"excitatory\"\n\
  14.     Gaussian - real positive number.\n\
  15. masksize [7] is the size of the mask (linear).\n\
  16. ratio [1.4] is the ratio b/w the s.d.'s of the inhibitory and\n\
  17.     excitatory Gaussians.\n\
  18. -p #    followed by a positive integer specifies the precision.\n\
  19.         Defaults to 1.\n\
  20. -b #    begin from #th frame. default is 0.\n\
  21. -n #    number of frames will be processed, default is total.\n\
  22. -s #    span -- frames between 2 slices.\n\
  23. -i    implies output in integer format.\n\
  24. -c    if -i is specified, causes checking of input to be in the\n\
  25.         range [-1024 to 1024].\n\
  26. -m    output the Gaussian(s) only.\n\
  27. -l    force dog3d do level x-y 2d filter    \n\
  28.         ([x is a level axle torward face,\n\
  29.         [ y axle is from left to right,    \n\
  30.         [ z is a vertical axle.])    \n\
  31. -v    force dog3d do vertical x-z 2d filter.    \n\
  32. -f    force dog3d do vertical y-z 2d filter.    \n\
  33. -g #    output the Gaussian(s) in filter_size + #\n\
  34.         Otherwise, output the Gaussian(s) in method ( 2**n ).\n\
  35. -w    output without header.\n\n\
  36.     If input file is not applied, the filtering is done on\n\
  37.         an impulse response in a 7 x 7 x 7 frame.\n";
  38. /*
  39. % compile:    cc -O -o $(DestDIR)/3dog 3dog.c -lccs -lhipsh -lhips -lm
  40. %        The -O Optimize option will speed up process.
  41. %
  42. % Author:    Guojun Jin - LBL    11/25/90
  43. */
  44.  
  45. #include "header.def"
  46. #include "imagedef.h"
  47. #include <math.h>
  48.  
  49. int    span;        /* frames in two slices    */
  50. MType    fsize;
  51. U_IMAGE    uimg;
  52.  
  53. #define    bbuf    uimg.src
  54. #define    col    uimg.width
  55. #define    row    uimg.height
  56. #define    frm    uimg.frames
  57. #define    in_fmt    uimg.in_form
  58.  
  59. #define    Float    float
  60. #define    Range    1024
  61. #define    Faster
  62.  
  63. #define    bgf    bgn_frm
  64. #define    k    num_offp
  65.  
  66. #ifdef    _DEBUG_
  67. extern    int    debug;
  68. #endif    _DEBUG_
  69.  
  70. #define    GValue()    arget(argc, argv, &i, &c)
  71.  
  72.  
  73. main(argc, argv)
  74. int    argc;
  75. char**    argv;
  76. {
  77. MType    *int_inbp, *int_obp, *int_conv,
  78.     *int_emask, *int_imask;
  79. int    synthetic,
  80.     woh=0, bgn_frm=1, num_offp=0,
  81.     isint=0, checking=0, gaussonly=0,
  82.     precision=1, masksize=7,
  83.     i, edglen,
  84.     ovf, unf;
  85. bool    Dir=0;
  86.  
  87.     /*    using same variable at different parse    */
  88. #define    form_size    gaussonly
  89. #define    c    edglen
  90.  
  91. byte    *bp;    /* byte buffer & pointer    */
  92. char    *fname=NULL;
  93. Float    *inbuf, *obuf,    /* buffers    */
  94.     *conv,        /* temp converting buffer    */
  95.     *obp, *inbp,    /* in & out put buffer pointers    */
  96.     *emask, *imask, gaussize=0;
  97. double    esigma=1., ratio=1.4,
  98.     gauss_mask();
  99.  
  100. format_init(&uimg, IMAGE_INIT_TYPE, HIPS, -1, *argv, "S20-1");
  101.  
  102. for (i = 1; i < argc; i++)
  103.   if (*argv[i] == '-') {
  104.     c = 1;
  105. args:    switch(argv[i][c++])    {
  106.     case 'p':
  107.         precision = GValue();    break;
  108.     case 'A':
  109.         esigma    = GValue();
  110.         if (!esigma)    goto    args;
  111.         i++;
  112.         frm = masksize= GValue();
  113.         if (!frm)    goto    args;
  114.         i++;
  115.         ratio    = arget(argc, argv,&i,&c,0);    break;
  116.     case 'i':
  117.         isint++;    break;
  118.     case 'c':
  119.         checking++;    break;
  120.     case 'b':
  121.         bgn_frm = GValue();    break;
  122.     case 'n':
  123.         num_offp = GValue();    break;
  124. #ifdef    _DEBUG_
  125.     case 'd':
  126.         debug = TRUE;    break;
  127. #endif    _DEBUG_
  128.     case 'm':
  129.         gaussonly++;
  130.         msg("%s no convolution; output the Gaussian(s)\n", Progname);
  131.         break;
  132.     case 's':
  133.         span = GValue() + 1;    break;
  134.     case 'l':
  135.     case 'v':
  136.     case 'f':
  137.         Dir = *argv[i];    break;
  138.     case 'g':
  139.         gaussize = GValue();    break;
  140.     case 'w':
  141.         woh++;    break;
  142.     default:
  143. info:        usage_n_options(usage, i, argv[i]);
  144.     }
  145.    }
  146.    else    fname = argv[i];
  147.  
  148. io_test(fileno(out_fp), goto    info);
  149.  
  150. if (esigma <= 0.0)    syserr("sigma must be positive");
  151. if (masksize < 1)    syserr("mask-size must be > 1");
  152. if (ratio < 0.0)    syserr("ratio must be >= 0");
  153. if (ratio == 1)        ratio = 0.0;
  154.  
  155. if (checking && !isint) {
  156.     mesg("checking mode ignored; output is not INT.\n");
  157.     checking=0;
  158. }
  159.  
  160. if (fname){
  161.     if ((in_fp=freopen(fname, "rb", stdin))==NULL)
  162.         syserr("input file %s", fname);
  163.     goto    rdh;
  164. }
  165. else if (!system("test -t 0")) {    /* no input file */
  166.     edglen = 1;
  167.     if (gaussize)
  168.         edglen = masksize + gaussize;
  169.     else    while(edglen < masksize)    edglen <<= 1;
  170.     if (frm > masksize)    frm = masksize>>1;
  171.     else    if (!frm)    frm = masksize;
  172.     if (!(frm % 2))    frm--;    /* must odd frames    */
  173.     if (Dir || frm<1)    frm = 1;
  174.     if (!gaussonly)
  175.         message("%s: output %dx%dx%d sample\n", *argv,edglen,edglen,frm);
  176.     col = row = edglen;
  177.     uimg.pxl_out = sizeof(float);
  178.     uimg.o_form = IFMT_FLOAT;
  179.     synthetic = 1;
  180.     (*uimg.std_swif)(FI_INIT_NAME, &uimg, "DOG mask", 0);
  181.     }
  182. else    {
  183. rdh:    (*uimg.header_handle)(HEADER_READ, &uimg, 0, 0);
  184.     if (uimg.in_form != IFMT_BYTE && uimg.in_form !=IFMT_FLOAT)
  185.         syserr("image pixel format must be bytes or floats");
  186.     if (uimg.frames < 2)
  187.         syserr("Can not handle non3D images");
  188.     synthetic = 0;
  189.     }
  190.  
  191. fsize = row * col;    /* size of one frame    */
  192. if (!bgn_frm || bgn_frm > frm)    bgn_frm=0;
  193. else    bgn_frm--;
  194. if (!num_offp || num_offp+bgn_frm > frm) num_offp = frm-bgn_frm;
  195. frm = num_offp;
  196.  
  197. if (!ratio)    msg("%s blurring by a Gaussian mask\n", Progname);
  198.  
  199. if (!gaussonly)
  200.    {
  201.     if (in_fmt==IFMT_BYTE)
  202.         bbuf = (VType *) zalloc(fsize, (MType)frm, "bbuf");
  203.     inbuf = zalloc(fsize, (MType)sizeof(Float) * frm, "inbuf");
  204.     conv = zalloc(fsize, (MType)sizeof(Float) * frm, "conv");
  205.     obuf = zalloc(fsize, (MType)sizeof(Float) * frm, "obuf");
  206.     uimg.o_form = IFMT_FLOAT;
  207.     uimg.pxl_out = sizeof(float);
  208.     if (isint) {
  209.         int_conv = (MType*)conv;
  210.         int_emask = zalloc((MType)masksize, (MType)sizeof(MType));
  211.         int_imask = zalloc((MType)masksize, (MType)sizeof(MType));
  212.         uimg.o_form = IFMT_LONG;
  213.     }
  214.     if (woh)    /*    go onto screen    */
  215.         (*uimg.header_handle)(HEADER_FWRITE, &uimg, argc, argv, stderr);
  216.     else    (*uimg.header_handle)(HEADER_WRITE, &uimg, argc, argv, True);
  217.    }
  218.  
  219. emask = (Float *) zalloc((MType)masksize, (MType)sizeof(Float), "emask");
  220. imask = (Float *) zalloc((MType)masksize, (MType)sizeof(Float), "imask");
  221.     gauss_mask(esigma,masksize,emask,precision,gaussonly?out_fp:NULL);
  222. if (ratio)
  223.     gauss_mask(esigma*ratio,masksize,imask,precision,gaussonly?out_fp:NULL);
  224. if (gaussonly)    return(0);    /* Done    */
  225.  
  226. if (isint)    {
  227.     for(i=0; i<masksize; i++) {
  228.         int_emask[i] = emask[i] * Range + .5;
  229.         int_imask[i] = imask[i] * Range + .5;
  230.     }
  231.     form_size = sizeof(MType);
  232. }
  233. else    form_size = sizeof(Float);
  234.  
  235. if (bgn_frm)    fseek(in_fp, (bgn_frm)*fsize, 1);
  236.  
  237. if (synthetic){
  238.     if (edglen & 1)    /* is an odd# ?    */
  239.         synthetic = (frm >> 1)*fsize + (fsize >> 1) + 1;
  240.     else    synthetic = (frm * fsize + col) >> 1;
  241.     if (isint) {
  242.         int_inbp = (MType *)inbuf;
  243.         for(i=0; i < fsize*frm; i++)    int_inbp[i] = 0;
  244.         int_inbp[synthetic] = 1; /* center=1 */
  245.     }
  246.     else {
  247.         for(i=0; i<fsize; i++)    inbuf[i] = 0.0;
  248.         inbuf[synthetic] = 1.0;
  249.     }
  250. }
  251. else    /* read whole image --> memory    */
  252.     if (in_fmt==IFMT_BYTE)
  253.     {    (*uimg.std_swif)(FI_LOAD_FILE, &uimg, uimg.load_all=frm, No);
  254.         if (isint)    /* type convert    */
  255.            for(i=0, bp=bbuf, int_inbp=(MType*)inbuf; i<fsize*frm; i++)
  256.             *int_inbp++ = *bp++ & 0xFF;
  257.         else for(i=0, bp=bbuf, inbp=inbuf; i<fsize*frm; i++)
  258.             *inbp++ = *bp++ & 0xFF;
  259.     }
  260.     else{
  261.         if (upread(inbuf, form_size, fsize * frm, in_fp) !=
  262.             fsize * frm)
  263. rerror:            syserr("unexpected EOF");
  264.         if (isint)     /* type convert; Float -> int */
  265.            for(i=0,int_inbp=(MType*)(inbp=inbuf); i<fsize*frm;i++)
  266.             *int_inbp++ = *inbp++;
  267.         }
  268.  
  269. if (checking)
  270.     for(i=unf=ovf=0,int_inbp=(MType*)inbuf; i<fsize*frm; i++){
  271.     k = *int_inbp++;
  272.     unf += (k < -Range);
  273.     ovf += (k > Range);
  274.     }
  275.     /*    starting convolution    */
  276.  
  277. if (isint) {
  278.     msg("%s_S4-1: integer (%c) GAUSSIAN\n", Progname, Dir);
  279.    if (!Dir)
  280.     {
  281.     int_vconvolve(int_emask,masksize,inbuf,obuf,row,col,frm);
  282.     int_hconvolve(int_emask,masksize,obuf,int_conv,row,col,frm);
  283.     int_fconvolve(int_emask,masksize,int_conv,obuf,row,col,frm);
  284.     if (ratio)
  285.        {    int_vconvolve(int_imask,masksize,inbuf,int_conv,row,col,frm);
  286.         int_hconvolve(int_imask,masksize,int_conv,inbuf,row,col,frm);
  287.         int_fconvolve(int_imask,masksize,inbuf,int_conv,row,col,frm);
  288.         free(inbuf);    inbuf = (Float*)int_conv;
  289.        }
  290.     }
  291.    else    switch(Dir)
  292.     {
  293.     case 'f':
  294.     int_vconvolve(int_emask,masksize,inbuf,int_conv,row,col,frm);
  295.     int_hconvolve(int_emask,masksize,int_conv,obuf,row,col,frm);
  296.     if (ratio)
  297.        {    int_vconvolve(int_imask,masksize,inbuf,int_conv,row,col,frm);
  298.         int_hconvolve(int_imask,masksize,int_conv,inbuf,row,col,frm);
  299.        }
  300.     break;
  301.     case 'l':
  302.     int_hconvolve(int_emask,masksize,inbuf,int_conv,row,col,frm);
  303.     int_fconvolve(int_emask,masksize,int_conv,obuf,row,col,bgf,frm);
  304.     if (ratio)
  305.        {    int_hconvolve(int_imask,masksize,inbuf,int_conv,row,col,frm);
  306.         int_fconvolve(int_imask,masksize,int_conv,obuf,row,col,frm);
  307.        }
  308.     break;
  309.     case 'v':
  310.     int_vconvolve(int_emask,masksize,inbuf,int_conv,row,col,frm);
  311.     int_fconvolve(int_emask,masksize,int_conv,obuf,row,col,frm);
  312.     if (ratio)
  313.        {    int_vconvolve(int_imask,masksize,inbuf,int_conv,row,col,frm);
  314.         int_fconvolve(int_imask,masksize,int_conv,obuf,row,col,frm);
  315.        }
  316.     break;
  317.     }
  318.    if (ratio)
  319.     for(i=0, int_obp=(MType*)obuf, int_inbp=(MType*)inbuf; i<fsize*frm; i++)
  320.         *int_obp++ -= *int_inbp++;
  321. }
  322. else{    msg("%s_S4-1: regular [%c span=%d] GAUSSIAN\n", Progname, Dir, span);
  323.    if (!Dir){
  324.     vconvolve(emask,masksize,inbuf,obuf,row,col,frm);
  325.     hconvolve(emask,masksize,obuf,conv,row,col,frm);
  326.     fconvolve(emask,masksize,conv,obuf,row,col,frm);
  327.     if (ratio) {
  328.         vconvolve(imask,masksize,inbuf,conv,row,col,frm);
  329.         hconvolve(imask,masksize,conv,inbuf,row,col,frm);
  330.         fconvolve(imask,masksize,inbuf,conv,row,col,frm);
  331.         inbuf = conv;
  332.     }
  333.    }
  334.    else    switch(Dir)
  335.     {
  336.     case 'f':
  337.         vconvolve(emask,masksize,inbuf,conv,row,col,frm);
  338.         hconvolve(emask,masksize,conv, obuf,row,col,frm);
  339.         if (ratio) {
  340.             vconvolve(imask,masksize,inbuf,conv,row,col,frm);
  341.             hconvolve(imask,masksize,conv,inbuf,row,col,frm);
  342.         }
  343.         break;
  344.     case 'l':
  345.         hconvolve(emask,masksize,inbuf,conv,row,col,frm);
  346.         fconvolve(emask,masksize,conv, obuf,row,col,frm);
  347.         if (ratio) {
  348.             hconvolve(imask,masksize,inbuf,conv,row,col,frm);
  349.             fconvolve(imask,masksize,conv,inbuf,row,col,frm);
  350.         }
  351.         break;
  352.     case 'v':
  353.         vconvolve(emask,masksize,inbuf,conv,row,col,frm);
  354.         fconvolve(emask,masksize,conv, obuf,row,col,frm);
  355.         if (ratio) {
  356.             vconvolve(imask,masksize,inbuf,conv,row,col,frm);
  357.             fconvolve(imask,masksize,conv,inbuf,row,col,frm);
  358.         }
  359.         break;
  360.     }
  361.    if (ratio)
  362.     for(i=0, obp=obuf, inbp=inbuf; i < fsize*frm; i++)
  363.         *obp++ -= *inbp++;
  364. }    /* end of convolution    */
  365. fwrite(obuf, frm * form_size, fsize, out_fp);
  366.  
  367. if (checking && (ovf || unf))
  368.     message("%s: %d overflows and %d underflows !!\n", Progname,ovf,unf);
  369. }
  370.  
  371. #ifdef    Faster
  372.  
  373. int_hconvolve(mask,masksize, inb, outb,r,c,f)
  374. MType    *mask, *inb, *outb;
  375. {
  376. MType    row_n,        /* row # -- position    */
  377.     v_lft_cln,    /* very left column    */
  378.     i, hrzn_p,    /* working point    */
  379.     *tmpp, *p, pc,    /* relative column pos    */
  380.     fn, *obp=outb,
  381.     vlpix, vrpix,
  382.     sum, lweight, rweight,
  383.  
  384. last_cln= c - 1,
  385. msk2    = masksize >> 1,
  386.  
  387. lftregion = msk2 > c ? c : msk2,
  388. mid_rgn = c - (masksize-msk2);    if (mid_rgn < lftregion) mid_rgn = 0;
  389.  
  390. for (fn=0; fn<f; fn++, inb+=fsize)
  391.   for (row_n=v_lft_cln=0; row_n < r; row_n++, v_lft_cln+=c)
  392.   {    tmpp = inb + v_lft_cln;    /* current row & first column    */
  393.     vlpix = *tmpp;        /* cur row very left side pixel    */
  394.     vrpix = *(tmpp + last_cln);    /* very rlght side pix    */
  395.     for(hrzn_p=0; hrzn_p < lftregion; hrzn_p++) {
  396.         pc = hrzn_p - msk2;
  397.         p = tmpp + pc;
  398.         /* PC moved form left half mask to right half mask    */
  399.         for(i=sum=lweight=rweight=0; i < masksize; i++,pc++,p++)
  400.            {
  401.             if (pc < 0)
  402.                 lweight += mask[i];
  403.             else    if (pc < c)
  404.                     sum += mask[i] * *p;
  405.             else    rweight += mask[i];
  406.            }
  407.         *obp++ = sum + vlpix*lweight + vrpix*rweight;
  408.     }
  409.     for(; hrzn_p <= mid_rgn; hrzn_p++) {
  410.         p = tmpp + hrzn_p - msk2;
  411.         for(i=sum=0; i < masksize;i++,p++)
  412.             sum += *p * mask[i];
  413.         *obp++ = sum;
  414.     }
  415.     for(; hrzn_p < c; hrzn_p++) {
  416.         pc = hrzn_p - msk2;
  417.         p = tmpp + pc;
  418.         for(i=sum=rweight=0; i < masksize; i++,pc++,p++) {
  419.             if (pc > last_cln)
  420.                 rweight += mask[i];
  421.             else    sum += *p * mask[i];
  422.         }
  423.         *obp++ = sum + vrpix*rweight;
  424.     }
  425.    }
  426. }
  427.  
  428. int_vconvolve(mask, masksize, inb, outb, r, c, f)
  429. MType    *mask, *inb, *outb;
  430. {
  431. MType    row_n, hrzn_p, i, fn,
  432.     sum,
  433.     bweight, tweight,
  434.     *tmpp, *p, pr, *obp = outb,
  435.  
  436. msk2    = masksize >> 1,
  437. last_row= r - 1,
  438. nd_r_p    = last_row * c,    /* last row (end row) first column    */
  439.  
  440. topregion = msk2 > r ? r : msk2,
  441. mid_rgn = r-(masksize-msk2); if (mid_rgn < topregion) mid_rgn = 0;
  442.  
  443. for (fn=0; fn < f; fn++, inb+=fsize) {
  444.   for (row_n=0; row_n < topregion; row_n++)
  445.     for (hrzn_p=0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++)
  446.     {
  447.     p = tmpp + hrzn_p;    pr = row_n - msk2;
  448.     for(i=sum=bweight=tweight=0; i<masksize; i++,p+=c,pr++)
  449.        {
  450.         if (pr < 0)        bweight += mask[i];
  451.         else    if (pr < r)    sum += *p * mask[i];
  452.             else    tweight += mask[i];
  453.        }
  454.     *obp++ = sum + inb[hrzn_p]*bweight + inb[hrzn_p + nd_r_p]*tweight;
  455.     }
  456.   for(; row_n <= mid_rgn; row_n++)
  457.     for(hrzn_p=0, tmpp=inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++) {
  458.     for(i=sum=0, p=tmpp + hrzn_p; i < masksize; i++, p+=c)
  459.         sum += *p * mask[i];
  460.     *obp++ = sum;
  461.     }
  462.   for(; row_n < r; row_n++)
  463.     for(hrzn_p=0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++) {
  464.     p = tmpp + hrzn_p;
  465.     for(i=sum=tweight=0,pr=row_n - msk2; i < masksize;i++,pr++,p+=c)
  466.         if (pr > last_row)
  467.             tweight += mask[i];
  468.         else     sum += *p * mask[i];
  469.     *obp++ = sum + tweight * inb[nd_r_p + hrzn_p];
  470.     }
  471.   }
  472. }
  473.  
  474. int_fconvolve(mask,masksize,inb,outb,r,c,f)
  475. MType    *mask, *inb,*outb;
  476. {
  477. MType    row_n, hrzn_p,
  478.     frm_n,
  479.     msk2 = masksize >> 1,
  480.     last_frm = f - 1,
  481.     i, pf, *tmpp,*p,/* current point (temp)    */
  482.     sum, weight,    /* rear */
  483.     fweight,    /* front */
  484.     *obp = outb,
  485.  
  486. frtregion = msk2 > f ? f : msk2,
  487. mid_rgn = f - (masksize-msk2);    if (mid_rgn < frtregion) mid_rgn = 0;
  488.  
  489. for (frm_n=0; frm_n < frtregion; frm_n++)
  490.   for (row_n=0, tmpp=inb + (frm_n - msk2)*fsize; row_n < r; row_n++,tmpp+=c)
  491.     for (hrzn_p=0; hrzn_p < c; hrzn_p++) {
  492.     p = tmpp + hrzn_p;
  493.     sum = fweight = weight = 0;
  494.     for(i=0; i<masksize; i++, p+=fsize) {
  495.         pf = frm_n - msk2 + i;
  496.         if (pf < 0)        fweight += mask[i];
  497.         else    if (pf < f)    sum += *p * mask[i];
  498.         else    weight += mask[i];
  499.     }
  500.     *obp++ = sum + inb[row_n*c+hrzn_p]*fweight +
  501.         inb[hrzn_p + (last_frm*r+row_n)*c]*weight;
  502.    }
  503. for (; frm_n < mid_rgn; frm_n++)
  504.   for (row_n = 0, tmpp=inb + (frm_n - msk2)*fsize; row_n < r; row_n++,tmpp+=c)
  505.     for (hrzn_p = 0; hrzn_p < c; hrzn_p++) {
  506.     p = tmpp + hrzn_p;
  507.     for(i = sum = 0; i < masksize; i++, p+=fsize)
  508.         sum += *p * mask[i];
  509.     *obp++ = sum;
  510.    }
  511. for (; frm_n < f; frm_n++)
  512.   for (row_n = 0, tmpp=inb + (frm_n - msk2)*fsize; row_n < r; row_n++,tmpp+=c)
  513.     for (hrzn_p = 0; hrzn_p < c; hrzn_p++) {
  514.     p = tmpp + hrzn_p;
  515.     sum = weight = 0;
  516.     for(i = 0; i < masksize; i++, p+=fsize) {
  517.         pf = frm_n - msk2 + i;
  518.         if (pf < f)    sum += *p * mask[i];
  519.         else    weight += mask[i];
  520.     }
  521.     *obp++ = sum + inb[hrzn_p + (last_frm*r+row_n)*c]*weight;
  522.    }
  523. }
  524.  
  525. hconvolve(mask,masksize, inb, outb,r,c,f)
  526. Float    *mask, *inb, *outb;
  527. {
  528. MType    row_n, hrzn_p, i, pc,
  529.     v_lft_cln, fn;
  530. Float    vlpix, vrpix,
  531.     sum, lweight, rweight,
  532.     *tmpp, *p,
  533.     *obp = outb;
  534.  
  535. int
  536. last_cln    = c - 1,
  537. msk2    = masksize >> 1,
  538.  
  539. lftregion= msk2 > c ? c : msk2,
  540. mid_rgn= c-(masksize-msk2); if (mid_rgn < lftregion) mid_rgn = 0;
  541.  
  542. for (fn = 0; fn < f; fn++, inb += fsize)
  543.   for (row_n = v_lft_cln = 0; row_n < r; row_n++, v_lft_cln += c)
  544.   {    tmpp = inb + v_lft_cln;    /* current row & first column    */
  545.     vlpix = *tmpp;    /* current row & very left side pixel    */
  546.     vrpix = *(tmpp + last_cln);    /* very rlght side pix    */
  547.     for(hrzn_p = 0; hrzn_p < lftregion; hrzn_p++) {
  548.         pc = hrzn_p - msk2;
  549.         p = tmpp + pc;
  550.         /* PC moved form left half mask to right half mask    */
  551.         for(i=sum=lweight=rweight = 0; i < masksize; i++,pc++,p++)
  552.            {
  553.             if (pc < 0)
  554.                 lweight += mask[i];
  555.             else    if (pc < c)
  556.                     sum += mask[i] * *p;
  557.             else    rweight += mask[i];
  558.            }
  559.         *obp++ = sum + vlpix*lweight + vrpix*rweight;
  560.     }
  561.     for(; hrzn_p <= mid_rgn; hrzn_p++) {
  562.         p = tmpp + hrzn_p - msk2;
  563.         for(i = sum = 0; i < masksize;i++,p++)
  564.             sum += *p * mask[i];
  565.         *obp++ = sum;
  566.     }
  567.     for(; hrzn_p < c; hrzn_p++) {
  568.         pc = hrzn_p - msk2;
  569.         p = tmpp + pc;
  570.         for(i = sum = rweight = 0; i < masksize; i++,pc++,p++) {
  571.             if (pc > last_cln)
  572.                 rweight += mask[i];
  573.             else    sum += *p * mask[i];
  574.         }
  575.         *obp++ = sum + vrpix*rweight;
  576.     }
  577.    }
  578. }
  579.  
  580. vconvolve(mask,masksize,inb,outb,r,c,f)
  581. Float    *mask, *inb,*outb;
  582. {
  583. MType    row_n, hrzn_p, i,
  584.     pr, fn;
  585. Float    sum , bweight , tweight,
  586.     *tmpp, *p,
  587.     *obp = outb;
  588.  
  589. MType    last_row = r - 1,
  590. msk2    = masksize >> 1,
  591. nd_r_p    = last_row * c,    /* last row (end row) first column    */
  592.  
  593. topregion= msk2 > r ? r : msk2,
  594. mid_rgn= r-(masksize-msk2); if (mid_rgn < topregion) mid_rgn = 0;
  595.  
  596. for (fn = 0; fn < f; fn++, inb += fsize) {
  597.   for (row_n = 0; row_n < topregion; row_n++)
  598.     for (hrzn_p = 0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++)
  599.     {
  600.     p = tmpp + hrzn_p;    pr = row_n - msk2;
  601.     for(i=sum=bweight=tweight=0; i<masksize; i++,p+=c,pr++)
  602.        {
  603.         if (pr < 0)        bweight += mask[i];
  604.         else    if (pr < r)    sum += *p * mask[i];
  605.             else    tweight += mask[i];
  606.        }
  607.     *obp++ = sum + inb[hrzn_p]*bweight + inb[hrzn_p + nd_r_p]*tweight;
  608.     }
  609.   for(; row_n <= mid_rgn; row_n++)
  610.     for(hrzn_p = 0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++) {
  611.     for(i=sum= 0, p = tmpp + hrzn_p; i < masksize; i++, p+=c)
  612.         sum += *p * mask[i];
  613.     *obp++ = sum;
  614.     }
  615.   for(; row_n < r; row_n++)
  616.     for(hrzn_p = 0, tmpp = inb + (row_n-msk2)*c; hrzn_p < c; hrzn_p++) {
  617.     p = tmpp + hrzn_p;
  618.     for(i=sum=tweight= 0,pr = row_n - msk2; i < masksize;i++,pr++,p+=c)
  619.         if (pr > last_row)
  620.             tweight += mask[i];
  621.         else     sum += *p * mask[i];
  622.     *obp++ = sum + tweight * inb[nd_r_p + hrzn_p];
  623.     }
  624.   }
  625. }
  626.  
  627. /*    interpolate from front to rear --
  628. *    when the frame interpolated is front of certain frame, start at
  629. *    more front one frame. If the frame is interpolated torwrad rear,
  630. *    then just start at that frame.
  631. */
  632. Float    frm_interpolation(masksize, mask, pp, span, fn, last_frm, Fv, Bv)
  633. MType    fn, last_frm;
  634. register Float    *mask, *pp;
  635. {
  636. MType    i, emf, realf, pf, num, Fw=0, Bw=0,
  637.     sfs = masksize / (span << 1),
  638.     smod = masksize % (span << 1);
  639. Float    d, val, sum=0;
  640.  
  641. if (smod-1)
  642.     sfs = (sfs + 1) << 1;
  643. else    sfs = (sfs << 1) + 1;
  644.  
  645. for (i=0; i < sfs; i++)
  646.   {
  647.     if (i && i < sfs-1)    num = span;
  648.     else if (!i)    num = smod >> 1;
  649.     else    num = smod - (smod >> 1);
  650.     if (!num)    num = span;
  651.  
  652.     pf =  i - (sfs >> 1);
  653.     realf = pf + fn;
  654.     if (realf<0)    for (emf=0; emf<num; emf++)    Fw += *mask++;
  655.     else if (realf >= last_frm)
  656.         for (emf=0; emf<num; emf++)
  657.         {    if (realf == last_frm && !emf)
  658.                 sum += *(pp + pf*fsize) * *mask++;
  659.             else    Bw += *mask++;
  660.         }
  661.     else {
  662.         val = *(pp + pf*fsize);
  663.         d = (*(pp + (pf + 1)*fsize) - val) / span;
  664.         for (emf=0; emf<num; emf++)
  665.             sum += (val + d*(span-num+emf)) * *mask++;
  666.     }
  667.   }
  668. return    Fw *  Fv + Bw * Bv + sum;
  669. }
  670.  
  671. fconvolve(mask,masksize,inb,outb,r,c,f)
  672. Float    *mask, *inb,*outb;
  673. {
  674. MType    row_n, hrzn_p,
  675.     frm_n,
  676.     msk2 = masksize >> 1,
  677.     last_frm = f - 1,
  678.     i, pf;
  679. Float    sum, weight,    /* rear */
  680.     fweight,    /* front */
  681.     *tmpp, *p,
  682.     *obp = outb;
  683.  
  684. MType    frtregion = msk2 > f ? f : msk2,
  685. mid_rgn = f - (masksize-msk2);    if (mid_rgn < frtregion) mid_rgn = 0;
  686.  
  687. if (span){
  688. int    Fv, Bv;
  689. /*
  690.     pf = (msk2 % span) != 0;
  691.     i = msk2 / span + pf;
  692.     frtregion = i > f ? f : i;
  693.     pf = ((masksize-msk2) % span) != 0;
  694.     mid_rgn = f - (masksize-msk2) / span + pf;
  695.     if (mid_rgn < frtregion) mid_rgn = 0;
  696.     else if (mid_rgn == f)    mid_rgn--;
  697.     wbuf = (Float*)zalloc((MType)masksize, (MType)sizeof(Float));
  698. */
  699.     msg("interpolate %d frames\n", span-1);
  700.  
  701. for (frm_n = 0; frm_n < f; frm_n++)
  702.   for (row_n = 0, tmpp=inb + frm_n*fsize; row_n < r; row_n++,tmpp+=c)
  703.     for (hrzn_p = 0; hrzn_p < c; hrzn_p++) {
  704.     p = tmpp + hrzn_p;
  705.     Fv = inb[row_n * c + hrzn_p];
  706.     Bv = inb[(last_frm*r + row_n) * c + hrzn_p];
  707.     *obp++ = frm_interpolation(masksize, mask, p, span, frm_n, last_frm, Fv, Bv);
  708.    }
  709. /****    faniest: if have tmpp = inb + (frm_n - msk2) * fsize    ****/
  710. /*
  711. for (; frm_n < mid_rgn; frm_n++)
  712.   for (row_n = 0, tmpp=inb + frm_n*fsize; row_n < r; row_n++,tmpp+=c)
  713.     for (hrzn_p = 0; hrzn_p < c; hrzn_p++) {
  714.     p = tmpp + hrzn_p;
  715.     if (frm_interpolation(masksize, wbuf, p, span, frm_n, last_frm))
  716.         syserr("interpolation frame %d", frm_n);
  717.     for(i = sum = 0; i < masksize; i++, p+=fsize)
  718.         sum += wbuf[i] * mask[i];
  719.     *obp++ = sum;
  720.    }
  721. for (; frm_n < f; frm_n++)
  722.   for (row_n=0, tmpp=inb + frm_n*fsize; row_n < r; row_n++,tmpp+=c)
  723.     for (hrzn_p=0; hrzn_p < c; hrzn_p++) {
  724.     p = tmpp + hrzn_p;
  725.     frm_interpolation(masksize, wbuf, p, span, frm_n, last_frm);
  726.     for(i=sum=weight=0; i < masksize; i++, p+=fsize)
  727.         if (wbuf > 0)    sum += wbuf[i] * mask[i];
  728.         else    weight += mask[i];
  729.     *obp++ = sum + inb[hrzn_p + (last_frm*r+row_n)*c]*weight;
  730.    }
  731. */
  732. }
  733. else    {
  734. for (frm_n=0; frm_n < frtregion; frm_n++)
  735.   for (row_n=0, tmpp=inb + (frm_n - msk2)*fsize; row_n < r; row_n++,tmpp+=c)
  736.     for (hrzn_p=0; hrzn_p < c; hrzn_p++) {
  737.     p = tmpp + hrzn_p;
  738.     sum = fweight = weight = 0;
  739.     for(i = 0; i < masksize; i++, p+=fsize) {
  740.         pf = frm_n - msk2 + i;    /* pf and p are consistent    */
  741.         if (pf < 0)        fweight += mask[i];
  742.         else    if (pf < f)    sum += *p * mask[i];
  743.         else    weight += mask[i];
  744.     }
  745.     *obp++ = sum + inb[row_n*c+hrzn_p]*fweight +
  746.         inb[hrzn_p + (last_frm*r+row_n)*c]*weight;
  747.    }
  748. for (; frm_n < mid_rgn; frm_n++)
  749.   for (row_n=0, tmpp=inb + (frm_n - msk2)*fsize; row_n < r; row_n++,tmpp+=c)
  750.     for (hrzn_p=0; hrzn_p < c; hrzn_p++) {
  751.     p = tmpp + hrzn_p;
  752.     for(i=sum=0; i < masksize; i++, p+=fsize)
  753.         sum += *p * mask[i];
  754.     *obp++ = sum;
  755.    }
  756. for (; frm_n < f; frm_n++)
  757.   for (row_n=0, tmpp=inb + (frm_n - msk2)*fsize; row_n < r; row_n++,tmpp+=c)
  758.     for (hrzn_p=0; hrzn_p < c; hrzn_p++) {
  759.     p = tmpp + hrzn_p;
  760.     for(i=sum=weight=0; i < masksize; i++, p+=fsize) {
  761.         pf = frm_n - msk2 + i;
  762.         if (pf < f)    sum += *p * mask[i];
  763.         else    weight += mask[i];
  764.     }
  765.     *obp++ = sum + inb[hrzn_p + (last_frm*r+row_n)*c]*weight;
  766.    }
  767. }
  768. }
  769.  
  770.  
  771. #else    /* easy to read and to understand    */
  772.  
  773. int_hconvolve(mask,masksize, inb, outb,r,c,f)
  774. MType    *mask, *inb, *outb;
  775. {
  776. MType    row_n,        /* row # -- position    */
  777.     v_lft_cln,    /* very left column    */
  778.     i, hrzn_p,    /* working point    */
  779.     msk2, pc,    /* relative column pos    */
  780.     fn, *obp=outb,
  781.     vlpix, vrpix,
  782.     sum, lweight, rweight;
  783.  
  784. msk2    = masksize >> 1;
  785.  
  786. for (fn = 0; fn < f; fn++, inb += fsize)
  787.   for (row_n=v_lft_cln=0; row_n < r; row_n++, v_lft_cln+=c)
  788.    {
  789.     vlpix = inb[v_lft_cln];        /* very left side pixel    */
  790.     vrpix = inb[v_lft_cln + c - 1];    /* very rlght side pix    */
  791.     for(hrzn_p=0; hrzn_p < c; hrzn_p++) {
  792.         pc = hrzn_p - msk2;
  793.         /*    PC moved form left half mask to right half mask    */
  794.         for(i = sum = lweight = rweight = 0; i < masksize; i++, pc++)
  795.            {
  796.             if (pc < 0)
  797.                 lweight += mask[i];
  798.             else    if (pc < c)
  799.                     sum += mask[i] * inb[v_lft_cln + pc];
  800.             else    rweight += mask[i];
  801.            }
  802.         *obp++ = sum + vlpix*lweight + vrpix*rweight;
  803.     }
  804.     }
  805. }
  806.  
  807. int_vconvolve(mask, masksize, inb, outb, r, c, f)
  808. MType    *mask, *inb, *outb;
  809. {
  810. MType    row_n, hrzn_p, i, fn,
  811.     sum , lweight, tweight,
  812.     pr, tp, *obp = outb,
  813.     msk2    = masksize >> 1,
  814. nd_r_p    = (r - 1) * c;    /* last row (end row) first column    */
  815.  
  816. for (fn = 0; fn < f; fn++, inb += fsize)
  817.   for (row_n = 0; row_n < r; row_n++)
  818.     for (hrzn_p = 0; hrzn_p < c; hrzn_p++)
  819.     {
  820.     for(i=sum=lweight=tweight=0, tp=(row_n-msk2)*c+hrzn_p;
  821.     i<masksize; i++,tp+=c)
  822.        {    pr = row_n - msk2 + i;
  823.         if (pr < 0)        lweight += mask[i];
  824.         else    if (pr < r)    sum += mask[i] * inb[tp];
  825.             else    tweight += mask[i];
  826.        }
  827.     *obp++ = sum + inb[hrzn_p]*lweight + inb[hrzn_p + nd_r_p]*tweight;
  828.     }
  829. }
  830.  
  831. int_fconvolve(mask,masksize,inb,outb,r,c,f)
  832. MType    *mask, *inb,*outb;
  833. {
  834. MType    row_n, hrzn_p,
  835.     frm_n,
  836.     msk2 = masksize >> 1,
  837.     last_frm = f - 1,
  838.     i, pf, tp,    /* current point (temp)    */
  839.     sum, weight,    /* rear */
  840.     fweight,    /* front */
  841.     *obp = outb;
  842.  
  843. for (frm_n=0; frm_n < f; frm_n++)
  844.   for (row_n=0; row_n < r; row_n++)
  845.     for (hrzn_p=0; hrzn_p < c; hrzn_p++) {
  846.     tp = ((frm_n - msk2)*r + row_n) * c + hrzn_p;
  847.     sum = fweight = weight = 0;
  848.     for(i = 0; i < masksize; i++, tp+=fsize) {
  849.         pf = frm_n - msk2 + i;
  850.         if (pf < 0)        fweight += mask[i];
  851.         else    if (pf < f)    sum += mask[i] * inb[tp];
  852.         else    weight += mask[i];
  853.     }
  854.     *obp++ = sum + inb[row_n*c+hrzn_p]*fweight +
  855.         inb[hrzn_p + (last_frm*r+row_n)*c]*weight;
  856.    }
  857. }
  858.  
  859. hconvolve(mask,masksize, inb, outb,r,c,f)
  860. Float    *mask, *inb, *outb;
  861. {
  862. MType    row_n, hrzn_p, i, pc,msk2,
  863.     v_lft_cln, fn;
  864. Float    vlpix, vrpix,
  865.     sum, lweight, rweight,
  866.     *obp = outb;
  867.  
  868. msk2    = masksize>>1;
  869.  
  870. for (fn=0; fn < f; fn++, inb+=fsize)
  871.   for (row_n=v_lft_cln=0; row_n < r; row_n++, v_lft_cln+=c)
  872.     {
  873.     vlpix = inb[v_lft_cln];
  874.     vrpix = inb[v_lft_cln + c - 1];
  875.     for(hrzn_p=0; hrzn_p < c; hrzn_p++) {
  876.         pc = hrzn_p - msk2;    /* start from over left to right */
  877.         for(i=sum=lweight=rweight=0; i < masksize; i++, pc++) {
  878.             if (pc < 0)    lweight += mask[i];
  879.             else    if (pc < c)
  880.                 sum += mask[i] * inb[v_lft_cln + pc];
  881.             else    rweight += mask[i];
  882.         }
  883.         *obp++ = sum + vlpix*lweight + vrpix*rweight;
  884.     }
  885.    }
  886. }
  887.  
  888. vconvolve(mask,masksize,inb,outb,r,c,f)
  889. Float    *mask, *inb,*outb;
  890. {
  891. MType    row_n, hrzn_p,i, pr, msk2,
  892.     nd_r_p, fn, tp;    /* current point (temp)    */
  893. Float    sum , bweight , tweight,
  894.     *obp = outb;
  895.  
  896. msk2    = masksize>>1;
  897. nd_r_p    = (r-1)*c;    /* last row (end row) first column    */
  898.  
  899. for (fn=0; fn < f; fn++, inb += fsize)
  900.   for (row_n=0; row_n < r; row_n++)
  901.     for (hrzn_p=0; hrzn_p < c; hrzn_p++) {
  902.     sum = bweight = tweight = 0;
  903.     for(i=0, tp = (row_n - msk2) * c + hrzn_p; i < masksize; i++, tp+=c)
  904.     {    pr = row_n - msk2 + i;
  905.         if (pr < 0)        tweight += mask[i];
  906.         else    if (pr < r)    sum += mask[i] * inb[tp];
  907.         else    bweight += mask[i];
  908.     }
  909.     *obp++ = sum + inb[hrzn_p]*tweight + inb[hrzn_p + nd_r_p]*bweight;
  910.    }
  911. }
  912.  
  913. fconvolve(mask,masksize,inb,outb,r,c,f)
  914. Float    *mask, *inb,*outb;
  915. {
  916. MType    row_n, hrzn_p,
  917.     frm_n,
  918.     msk2 = masksize>>1,
  919.     last_frm = f - 1,
  920.     i, pf, tp;    /* current point (temp)    */
  921. Float    sum, weight,    /* rear */
  922.     fweight,    /* front */
  923.     *obp = outb;
  924.  
  925. for (frm_n=0; frm_n < f; frm_n++)
  926.   for (row_n=0; row_n < r; row_n++)
  927.     for (hrzn_p=0; hrzn_p < c; hrzn_p++) {
  928.     tp = ((frm_n - msk2)*r + row_n) * c + hrzn_p;
  929.     sum = fweight = weight = 0;
  930.     for(i = 0; i < masksize; i++, tp+=fsize) {
  931.         pf = frm_n - msk2 + i;
  932.         if (pf < 0)        fweight += mask[i];
  933.         else    if (pf < f)    sum += mask[i] * inb[tp];
  934.         else    weight += mask[i];
  935.     }
  936.     *obp++ = sum + inb[row_n*c+hrzn_p]*fweight +
  937.         inb[hrzn_p + (last_frm*r+row_n)*c]*weight;
  938.    }
  939. }
  940. #endif    Faster
  941.